
/*=============================================================================================
																			
	Replication file for survey and price analysis:
	Towards Smallholder Food and Water Security

	William A. Sundstrom
	Current version: 12/4/2020

**=============================================================================================

	This file contains Stata code to replicate the following tables and figures 
	for the survey and price series tables and figures in the paper, 
	including supplementary material:
	
		Table 3: Key livelihood characteristics of matched survey households, 2014-2017
		Table 4: Drivers of food and water insecurity, 2017 survey: regressions 
		Figure 3: Relationship between measures of seasonal water insecurity and 
			food insecurity across survey sample households, 2017 (3-plot panel)
		Figure 4: Farm size and food insecurity across households, 2017
		Figure 5a: Intra-annual incidence of food and water insecurity, 
			household sample, 2017 (plot using R)
		Figure 5d: Seasonality of real basic grain prices, Nicaragua: 
			mean monthly deviation from linear trend (plot using R)
		Figure 6b: Ranking of adverse events in terms of impact during the past 
			12 months, 2017 household survey
		Figure 6c: Index of monthly prices of maize and beans relative to 
			price of coffee at annual harvest (January 2010 = 100) (plot using R)	
		Figure 6d: Export price of coffee (other milds), 1990-2018 (plot using R)	
		Appendix Table A1: Food and water insecurity coping responses, 2017 survey
		Appendix Tables A2 and A3: Alternative regression specifications
		Appendix Table A4: Drivers of income and crop production across households, 2017
		Additional calculations and alternative specifications mentioned in the text
	
	Required input data files (included in download):
		Survey: hh_survey_data.dta
		Price series: price_series.dta

	To run the code, follow the instructions to download and decompress the 
	root folder containing this file and the data files. 
	
	Edit line 61 below for the path to your root folder. 
			
*/

*** Settings 

	version 14
	clear all
	set more off
	set scheme s1color
	
	* ado files used: install if necessary
	* ssc install ietoolkit
	* net install estout, from("http://fmwww.bc.edu/RePEc/bocode/e") replace
	* net install st0563.pkg
	
*** Global macros to folders and data sets

	* Replace the following path with path to your root folder
	global root		"$root3" 
	global data		"$root/data"
	global survey	"$data/hh_survey_data.dta"
	global price 	"$data/price_series_data.dta"
	global output 	"$root/output"
	
    cd "$root/scrap" 

	set more off
	
	
*===============================================================================
*	Tables for main text of paper 
*===============================================================================
	
*** Set up the survey data

	set more off
	cap eststo clear
	
	* New variables for both 2014 and 2017
	use "$survey", clear
	
	* Set leaf rust impact to zero if missing (non-coffee producers)
	replace roya_impact_perc_num = 0 if roya_impact_perc_num==.
		
	* Rescale some variables for better tables
	* Precip in meters
	g ppt_m = ppt_14_mo/1000
	label var ppt_m "Community 14-mo total precip (meters)"
	* Fraction of coffee affected by roya
	g roya_impact_prop_num = roya_impact_perc_num/100
	label var roya_impact_prop_num "Proportion coffee affected by CLR"
	* crop production in 1000 kg
	g corn_1000kg = corn_kg_z / 1000
	g beans_1000kg = beans_red_kg_z /1000 
	g coffee_1000kg = coffee_kg_z / 1000
	label var corn_1000kg "Corn production last 12 mo. (1000 kg)"
	label var beans_1000kg "Bean production last 12 mo. (1000 kg)"
	label var coffee_1000kg "Coffee production last 12 mo. (1000 kg)"
	* Income in $1000
	g income_sum_dol_1000 = income_sum_dol/1000
	label var income_sum_dol_1000 "Household cash income in $1000"
	
	* Income per adult equivalent
	g income_sum_dol_pc_1000 = income_sum_dol_pc / 1000
	
	* log variables use asinh (inverse hyperbolic sin) to deal with small number of 0 values
	g log_income_sum_dol = asinh(income_sum_dol)
	label var log_income_sum_dol "Log (asinh) of household income"
	g log_farm_size_calc_ha = asinh(farm_size_calc_ha)
	label var log_farm_size_calc_ha "Log (asinh) of farm size"
	
	* Obtain 2014 educ characteristics for matched HHs (obs in both 2014 and 2017)
	
	* HH ID number
	sort hh_hash 
	egen hh_index = group(hh_hash)
	xtset hh_index time

	* assign 2014 education and respondent age to 2017 for matched HHs
	local key educ_years_head educ_sec_plus_head educ_years_max_hh respondent_age 
	foreach x in `key' { 
		replace `x' = l.`x' if year==2017
		} 
	replace respondent_age = respondent_age+3  if year==2017
	
	save temp_1417, replace

*===============================================================================
*	T3. Key livelihood characteristics of matched survey households, 2014-2017 
*===============================================================================

	set more off
	
	use temp_1417, clear
	
	g organic_producer = .
	replace organic_producer = coffee_organic if coffee_pos==1
	label var organic_producer "Produces organic coffee - missing if no coffee"
	
	local compare_short /// 
		num_hh_tot num_hh_under_15 farm_size_calc_ha food_grown_gt_half    /// 
		grows_basic_grains corn_kg_z corn_kg_p beans_red_kg_z beans_red_kg_p /// 
		trees_total_sum_z coffee_pos coffee_prod_or_dev /// 
		coffee_kg_z coffee_kg_p organic_producer roya_impact_gt_half_miss /// 
		income_sum_dol sells_coffee sells_corn sells_beans wage_sal_worker /// 
		food_lean_mo dietary_berry water_lean_mo ppt_m 
		
	* Use iebaltab to compare means across time
	* std means std instead of SE 
	iebaltab `compare_short', grpvar(year) rowvarlabels order(2017 2014) /// 
		std save("$output/T3_output.xlsx")  /// 
		tblnote("Compare means of selected variables 2014 vs 2017") replace 
		
	* obtain and save medians by year as well
	collapse (median) `compare_short', by(year)	
		
	export excel using "$output/T3_meds_output.xlsx", firstrow(var) replace 
	
*===============================================================================
*	T4. Drivers of food and water insecurity, 2017 survey: regressions 
*===============================================================================

*** Regressions 

	set more off
	cap eststo clear

	use temp_1417 if year==2017, clear
	
	* summary stats for dependent variables
	sum food_lean_mo food_coping_score water_any_lean_mo water_coping_score 
	sum food_coping_score if food_coping_score>0
	sum water_coping_score if water_coping_score>0
	sum farm_size_calc_ha ppt_m educ_sec_plus_head roya_impact_prop_num /// 
		income_sum_dol_1000 corn_1000kg beans_1000kg coffee_1000kg  
		
	* main regressors
	local regr farm_size_calc_ha educ_sec_plus_head ppt_m roya_impact_prop_num 
	
	*** food lean months
	* farm size, human K, precip, roya
	eststo, ti(food lean mo): reg food_lean_mo `regr', r
	* add income
	eststo, ti(food lean mo): reg food_lean_mo `regr' income_sum_dol_1000, r
	* replace income with crop production
	eststo, ti(food lean mo): reg food_lean_mo `regr' coffee_1000kg corn_1000kg beans_1000kg, r
	
	*** food coping
	* farm size, human K, precip, roya
	eststo, ti(food coping): tobit food_coping_score `regr', ll(0) vce(r)
	* add income
	eststo, ti(food coping): tobit food_coping_score `regr' income_sum_dol_1000, ll(0) vce(r)
	* replace income with crop production
	eststo, ti(food coping): tobit food_coping_score `regr' coffee_1000kg corn_1000kg beans_1000kg, ll(0) vce(r)

	*** water lean months
	* farm size, human K, precip, roya
	eststo, ti(water lean mo): reg water_any_lean_mo `regr', r
	* add income
	eststo, ti(water lean mo): reg water_any_lean_mo `regr' income_sum_dol_1000, r
	* replace income with crop production
	eststo, ti(water lean mo): reg water_any_lean_mo `regr' coffee_1000kg corn_1000kg beans_1000kg, r
		
	*** water coping
	* farm size, human K, precip, roya
	eststo, ti(water coping): tobit water_coping_score `regr', ll(0) vce(r)
	* add income
	eststo, ti(water coping): tobit water_coping_score `regr' income_sum_dol_1000, ll(0) vce(r)
	* replace income with crop production
	eststo, ti(water coping): tobit water_coping_score `regr' coffee_1000kg corn_1000kg beans_1000kg, ll(0) vce(r)
	
*** Write the table	
	esttab  using "$output/T4_output.csv", /// 
		se replace nonote nonumbers r2 pr2 mtitles nogaps label ///
		order(farm_size_calc_ha /// 
			educ_sec_plus_head ppt_m roya_impact_prop_num income_sum_dol_1000 /// 
			corn_1000kg beans_1000kg coffee_1000kg) /// 
		star(* .1 ** 0.05 *** 0.01)  /// 
		title("Food & water insecurity regressions")  ///
		addnotes("Notes: * p<0.1 ** p<0.05 *** p<0.01.") 
	eststo clear
	
	
*===============================================================================
*	Figures for main text of paper 
*===============================================================================
		
*===============================================================================
*	F3. Scatter plots: Relationship between water and food insecurity, 2017
*===============================================================================

*** Scatter plots for HH sample 2017 

	set more off
	use temp_1417 if year==2017, clear

	* Lean months including outliers
	* make the weights for the cell sizes
	bysort food_lean_mo water_any_lean_mo: egen food_water_cell_count = count(water_any_lean_mo)
	twoway (lfitci food_lean_mo water_any_lean_mo ) /// 
		(scatter food_lean_mo water_any_lean_mo [w=food_water_cell_count], /// 
		title("a. Lean months 2017", size(large)) /// 
		ytitle("Number of lean months – food", size(large)) /// 
		xtitle("Number of lean months – water", size(large) ) xscale(range(0 12)) xlabel(0(2)12) /// 
		xlabel(,grid labsize(large)) ylabel(,grid labsize(large)) /// 
		msymbol(circle_hollow) legend(off))	
	graph rename A, replace
	* graph export "$output/F3a_output.png", replace width(800) height(600)	
	sum food_lean_mo water_any_lean_mo if (food_lean_mo~=. & water_any_lean_mo~=.)
	reg food_lean_mo water_any_lean_mo, r
		
	* Coping scores
	* make the weights for the cell sizes
	bysort food_coping_score water_coping_score: egen coping_cell_count = count(water_coping_score)
	* Coping scores for HHs with positive scores on both 
	twoway (lfitci food_coping_score water_coping_score if food_coping_score>0 & water_coping_score>0) /// 
		(scatter food_coping_score water_coping_score [w=coping_cell_count] if food_coping_score>0 & water_coping_score>0, /// 
		title("b. Coping index 2017", size(large)) /// 
		ytitle("Coping score – food", size(large)) /// 
		xtitle("Coping score – water", size(large)) xscale(range(0 35)) xlabel(0(5)35) yscale(range(0 60)) ylabel(0(10)60)   /// 
		xlabel(,grid  labsize(large)) ylabel(,grid labsize(large)) ///  
		msymbol(circle_hollow) legend(off))	
	graph rename B, replace 
	* graph export "$output/F3b_output.png", replace width(800) height(600)	
	sum food_coping_score water_coping_score if /// 
		(food_coping_score~=. & water_coping_score~=. & food_coping_score>0 & water_coping_score>0)
	reg food_coping_score water_coping_score if food_coping_score>0 & water_coping_score>0, r
	
*** Scatter plot in changes
		
	* set up the cross-section in changes 

	set more off	
	use temp_1417, clear
	
	* Create changes in desired vars
	
	local key /// 
		food_lean_mo water_lean_mo   

	foreach x in `key' { 		
		g d_`x' = d.`x'
		_crcslbl d_`x' `x' 
		}
		
	* we now just keep 2017, having the d_values
	keep if time==1 
	keep d_food_lean_mo d_water_lean_mo 
	
	* Now the plot
	* make the weights for the cell sizes
	bysort d_food_lean_mo d_water_lean_mo : egen change_cell_count = count(d_food_lean_mo)
	
	* Lean months changes
	twoway (lfitci d_food_lean_mo d_water_lean_mo ) /// 
		(scatter d_food_lean_mo d_water_lean_mo [w=change_cell_count], /// 
		title("c. Change in lean months 2014-2017", size(large)) /// 
		ytitle("Change in lean months – food", size(large)) /// 
		xtitle("Change in lean months – water", size(large))  /// 
		xlabel(,grid labsize(large)) ylabel(,grid labsize(large))  /// 
		msymbol(circle_hollow) legend(off))	
	graph rename C, replace 
	* graph export "$output/F3c_output.png", replace width(800) height(600)	
	reg d_food_lean_mo d_water_lean_mo, r
	
***	Combined figure for paper
	
	graph combine A B C, cols(3) xsize(7.5) ysize(2.5)
	graph export "$output/F3panel_output.png", replace width(3600) height(1200)	
		

*===============================================================================
*	F4. Farm size and food insecurity across households, 2017
*===============================================================================
	
	use temp_1417 if year==2017, clear
	
	* Note farm size has a few very large outliers
	sum farm_size_calc_ha, d
	list farm_size_calc_ha if farm_size_calc_ha>=20	
	
	* remove farm size outliers >20 ha, scale circles to N
	bysort food_lean_mo farm_size_calc_ha: egen food_farm_cell_count = count(farm_size_calc_ha) if farm_size_calc_ha<=20
	twoway (lowess food_lean_mo farm_size_calc_ha if farm_size_calc_ha<=20, /// 
		lwidth(thick)) /// 
		(scatter food_lean_mo farm_size_calc_ha  [w=food_farm_cell_count] if farm_size_calc_ha<=20, /// 
		ytitle("Number of lean months – food")  xlabel(,grid) ylabel(,grid) /// 
		xtitle("Farm area (ha)" ) legend(off) xmtick(0(1)20, ticks) /// 
		msymbol(circle_hollow) mcolor(navy)) 
	graph export "$output/F4_output.png", replace width(800) height(600) 	
				
				
*===============================================================================
*	F5. Seasonality of food and water insecurity, precipitation, and agriculture 
*===============================================================================
		
*===============================================================================
*	F5a. Intra-annual incidence of food and water insecurity, 2017 
*===============================================================================

*** Plot seasonality of water and food shortage (by months)

	* food or water for drink cook clean or bathe
	
	set more off
	
	clear
	input mo str3 month
	1 "Jan"
	2 "Feb"
	3 "Mar"
	4 "Apr"
	5 "May"
	6 "Jun"
	7 "Jul"
	8 "Aug"
	9 "Sep"
	10 "Oct"
	11 "Nov"
	12 "Dec"
	end
	save months_temp, replace
	
	use temp_1417 if year==2017, clear
	
	foreach m in jan feb march april may june july aug sept oct nov dec { 
		g water_short_any_`m'_d = water_short_drink_`m'==1 | water_short_cook_`m'==1 | /// 
			water_short_clean_`m'==1 | water_short_bathe_`m'==1
		g lean_month_any_`m'_d = lean_month_`m'_d==1 | water_short_any_`m'_d==1
		local a lean_month_`m'_d water_short_any_`m'_d lean_month_any_`m'_d 
		cap local a_v `a_v' `a'
		}
	
	collapse `a_v' 
	save temp_a, replace

	set more off
	cap erase temp_mo.dta
	local mo = 0
	foreach m in jan feb march april may june july aug sept oct nov dec { 
		local mo = `mo' + 1
		use temp_a, clear
		g mo = `mo'
		g str month = "`m'"
		g p_short_food = lean_month_`m'_d 
		g p_short_water = water_short_any_`m'_d 
		g p_short_any = lean_month_any_`m'_d
		cap append using temp_mo
		save temp_mo, replace
		}

	* Save data for plotting in R
	keep mo p_short_food p_short_water p_short_any    
	merge 1:1 mo using months_temp
	drop _merge
	save "$output/F5a_data.dta", replace
		
	
*===============================================================================
*	F5d. Seasonality of basic grain prices, Nicaragua
*===============================================================================
			
*** Data for plotting seasonal pattern

	clear
	input mm str3 month
	1 "Jan"
	2 "Feb"
	3 "Mar"
	4 "Apr"
	5 "May"
	6 "Jun"
	7 "Jul"
	8 "Aug"
	9 "Sep"
	10 "Oct"
	11 "Nov"
	12 "Dec"
	end
	g mt = 683          /* fix prediction month to Dec 2016 */
	save mmm, replace
	
*** monthly seasonal effects for prices based on 2000-2016 averages

	* regressions with month dummies and linear time trend for Jan 2000-Dec 2016
	use "$price", clear	
	keep p_w_nica_beansred_real p_w_nica_maize_real mt mm yy /// 
		tt_maize_coffee_jf tt_beansred_coffee_jf 
		
	reg p_w_nica_beansred_real mt b13.mm if mt<684, noconstant         /* all years */
	use mmm, clear
	predict price
	g beans = price
	keep beans month mm
	save temp_bean_season, replace
	
	use "$price", clear		
	reg p_w_nica_maize_real mt b13.mm if mt<684, noconstant          /* all years */
	use mmm, clear
	predict price
	g maize = price
	keep maize month mm
	save temp_maize_season, replace
	
	use temp_bean_season, clear
	merge 1:1 mm using temp_maize_season
	drop _merge
	rename beans b
	rename maize m
	* Percentage difference relative to mean
	egen bbbb = mean(b) 
	egen mmmm = mean(m) 
	g beans = 100*b/bbbb-100
	g maize = 100*m/mmmm-100
		
	* Save data for plotting in R
	keep mm month maize beans 
	save "$output/F5d_data.dta", replace

	
*===============================================================================
*	F6b. Ranking of hazards in terms of impact during the past 12 months, 2017
*===============================================================================
	
	* Ranking of hazards in terms of impact during past 12 months
	
	set more off
	
	* count gives us total non-missing N
	use temp_1417 if year==2017, clear
	collapse (count) sever_* 

	* plain collapse is the mean for non-missing, so proportion naming at that rank	
	use "$survey" if year==2017, clear
	collapse sever_* 
	save temp_sever, replace
	
	* create the table	
	cap erase temp_tab4.dta
	foreach x in roya drought incr_price_food drop_price_crop incr_cost_agr flood /// 
			sickness loss_lvstock prop_loss death damg_home other {
		use temp_sever, clear
		keep sever_`x'_*
		rename sever_`x'_* sever_*
		g str hazard = "`x'"
		cap append using temp_tab4
		save temp_tab4, replace
		} 
		
	gsort -sever_top3
	order hazard sever_1 sever_2 sever_3 sever_top3
			
	export excel using "$output/F6b_output.xlsx", /// 
		firstrow(varlabels) sheet("raw") sheetreplace 
		
	
*===============================================================================
*	F6cd. Terms of trade between coffee and basic grains, 2000-2018
*===============================================================================

*** Time series price data

	use "$price", clear
	
	* Save data for plotting in R
	keep mt mm yy tt_maize_coffee_jf tt_beansred_coffee_jf OtherMilds
	save "$output/F6cd_data.dta", replace	
	
	
*===============================================================================
*	Appendix Tables
*===============================================================================
	
*===============================================================================
*	A1. Food and water insecurity coping responses, 2017
*===============================================================================

	set more off
		
	use temp_1417 if year==2017, clear
		 
	local food_coping_vars /// 	
		food_cope_dum_eat_dislike /// 	
		food_cope_dum_tortilla_supp /// 	
		food_cope_dum_eat_taro /// 	
		food_cope_dum_sell_fruit /// 	
		food_cope_dum_beans_last /// 	
		food_cope_dum_borrow_aid /// 	
		food_cope_dum_buy_credit /// 	
		food_cope_dum_wild_food /// 	
		food_cope_dum_sell_store_seed /// 		
		food_cope_dum_eat_elsewhere /// 	
		food_cope_dum_smaller_portions /// 	
		food_cope_dum_kids_eat /// 	
		food_cope_dum_less_meals /// 	
		food_cope_dum_skip_meals /// 	
		food_cope_dum_sell_resrc /// 	
		food_cope_dum_migrate 
		
	local water_coping_vars /// 
		water_cope_dum_bed_thirsty /// 
		water_cope_dum_no_all_day /// 
		water_cope_dum_not_wash_food /// 
		water_cope_dum_less_pref_srce /// 
		water_cope_dum_bathe_less /// 
		water_cope_dum_discuss /// 
		water_cope_dum_borrow /// 
		water_cope_dum_take /// 
		water_cope_dum_less_wash_cloth /// 
		water_cope_dum_less_cook_water /// 
		water_cope_dum_clean_home  /// 
		water_cope_dum_bathe_stream /// 
		water_cope_dum_less_wash_hands 
				
	collapse `food_coping_vars' `water_coping_vars'
	export excel using "$output/A1_output.xlsx", /// 
		firstrow(varlabels) sheet("raw") sheetreplace 
	
*===============================================================================
*	A2 & A3. Alternative regression specifications
*===============================================================================

	eststo clear

	* main regressors
	local regr farm_size_calc_ha educ_sec_plus_head ppt_m roya_impact_prop_num 
	
*** Poisson count regression of lean months

	use temp_1417 if year==2017, clear
	
	* food lean months
	eststo, ti(food poisson): poisson food_lean_mo `regr', vce(robust)
	eststo, ti(food poisson): poisson food_lean_mo `regr' income_sum_dol_1000, vce(robust)
	
	* water lean months
	eststo, ti(water poisson): poisson water_any_lean_mo `regr', vce(robust)
	eststo, ti(water poisson): poisson water_any_lean_mo `regr' income_sum_dol_1000, vce(robust)
		
*** Income with income per capita, log land area and log income per adult equivalent
		
	* food lean months income pc
	eststo, ti(food lean): reg food_lean_mo farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num income_sum_dol_pc_1000, r
	* food lean months log farm area and income 
	eststo, ti(food lean): reg food_lean_mo log_farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num log_income_sum_dol, r
	* food coping income pc
	eststo, ti(food coping): tobit food_coping_score farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num income_sum_dol_pc_1000, ll(0) vce(r)
	* food coping log farm area and income 
	eststo, ti(food coping): tobit food_coping_score log_farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num log_income_sum_dol, ll(0) vce(r)
		
	* food lean months income pc
	eststo, ti(water lean): reg water_any_lean_mo farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num income_sum_dol_pc_1000, r
	* food lean months log farm area and income 
	eststo, ti(water lean): reg water_any_lean_mo log_farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num log_income_sum_dol, r
	* food coping income pc
	eststo, ti(water coping): tobit water_coping_score farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num income_sum_dol_pc_1000, ll(0) vce(r)
	* food coping log farm area and income 
	eststo, ti(water coping): tobit water_coping_score log_farm_size_calc_ha educ_sec_plus_head /// 
		ppt_m roya_impact_prop_num log_income_sum_dol, ll(0) vce(r)

*** Table A2 	

	esttab  using "$output/A2_output.csv", /// 
		se replace nonote nonumbers r2 pr2 mtitles nogaps label ///
		order(farm_size_calc_ha log_farm_size_calc_ha /// 
			educ_sec_plus_head ppt_m roya_impact_prop_num /// 
			income_sum_dol_1000 income_sum_dol_pc_1000 log_income_sum_dol) /// 
		star(* .1 ** 0.05 *** 0.01)  /// 
		title("Food & water insecurity regressions")  ///
		addnotes("Notes: * p<0.1 ** p<0.05 *** p<0.01.") 
	eststo clear
	
*** Additional regressors

	local iregr farm_size_calc_ha educ_sec_plus_head ppt_m roya_impact_prop_num income_sum_dol_1000 
	local mregr farm_size_calc_ha educ_sec_plus_head ppt_m roya_impact_prop_num income_sum_dol_1000 /// 
		female respondent_age num_hh_tot num_hh_under_15   
	
	* add organizational affiliation dummies
	eststo, ti(food lean mo): reg food_lean_mo `iregr' org_1 org_2, r
	eststo, ti(food coping): tobit food_coping_score `iregr' org_1 org_2, ll(0) vce(r)
	eststo, ti(water lean mo): reg water_any_lean_mo `iregr' org_1 org_2, r
	eststo, ti(water coping): tobit water_coping_score `iregr' org_1 org_2, ll(0) vce(r)
	* add extended HH demographics	
	eststo, ti(food lean mo): reg food_lean_mo `mregr', r
	eststo, ti(food coping): tobit food_coping_score `mregr', ll(0) vce(r)
	eststo, ti(water lean mo): reg water_any_lean_mo `mregr', r
	eststo, ti(water coping): tobit water_coping_score `mregr', ll(0) vce(r)
			
*** Regressions in changes
		
	* set up the cross-section in changes 
	set more off	
	use temp_1417, clear
	* Create changes in desired vars
	local key /// 
		food_lean_mo water_lean_mo farm_size_calc_ha ppt_m /// 
		roya_impact_prop_num income_sum_dol_1000
	foreach x in `key' { 		
		g d_`x' = d.`x'
		_crcslbl d_`x' `x' 
		drop `x'
		rename d_`x' `x'
		}
	* we now just keep 2017, having the d_values
	keep if time==1 
	local iregr farm_size_calc_ha ppt_m roya_impact_prop_num income_sum_dol_1000 
	eststo, ti(ch food lean mo): reg food_lean_mo `iregr', r
	eststo, ti(ch water lean mo): reg water_lean_mo `iregr', r		
		
*** Table A3	

	esttab  using "$output/A3_output.csv", /// 
		se replace nonote nonumbers r2 pr2 mtitles nogaps label ///
		order(farm_size_calc_ha log_farm_size_calc_ha /// 
			educ_sec_plus_head ppt_m roya_impact_prop_num income_sum_dol_1000  /// 
			ft_org cac_org female respondent_age num_hh_tot num_hh_under_15) /// 
		star(* .1 ** 0.05 *** 0.01)  /// 
		title("Food & water insecurity regressions")  ///
		addnotes("Notes: * p<0.1 ** p<0.05 *** p<0.01.") 
	eststo clear

				
*===============================================================================
*	A4. Drivers of income and crop production across households, 2017
*===============================================================================
	
*** Production and income regressions

	set more off
	cap eststo clear

	use temp_1417 if year==2017, clear
		
	* income on farm size, human K, precip, roya
	eststo, ti(income): reg income_sum_dol_1000 farm_size_calc /// 
		educ_sec_plus_head ppt_m roya_impact_prop_num, r
	* coffee production	(missing = 0) on farm size, human K, precip, roya
	eststo, ti(coffee): reg coffee_kg_z farm_size_calc_ha /// 
		educ_sec_plus_head ppt_m roya_impact_prop_num, r
	* coffee production	(drop missing) on farm size, human K, precip, roya
	eststo, ti(coffee): reg coffee_kg_p farm_size_calc_ha /// 
		educ_sec_plus_head ppt_m roya_impact_prop_num, r
	* corn production	(missing = 0)
	eststo, ti(corn): reg corn_1000kg farm_size_calc_ha /// 
		educ_sec_plus_head ppt_m roya_impact_prop_num, r
	* beans production	(missing = 0)
	eststo, ti(beans): reg beans_1000kg farm_size_calc_ha /// 
		educ_sec_plus_head ppt_m roya_impact_prop_num, r
	* income on farm size, human K, precip, roya
	eststo, ti(income): reg income_sum_dol_1000 farm_size_calc /// 
		educ_sec_plus_head ppt_m roya_impact_prop_num, r
		
	*** Write the table	
	esttab  using "$output/A4_output.csv", /// 
		se replace nonote nonumbers r2 mtitles nogaps label ///
		order(farm_size_calc_ha /// 
			educ_sec_plus_head ppt_m roya_impact_prop_num) /// 
		star(* .1 ** 0.05 *** 0.01)  /// 
		title("Income and crop production regressions")  ///
		addnotes("Notes: * p<0.1 ** p<0.05 *** p<0.01.") 
	eststo clear

	
*===============================================================================
*	Additional calculations and specifications mentioned in the paper
*===============================================================================

*** Coffee transitions 2014-2017

	* set up the cross-section in changes 

	use "$survey", clear
	
	* HH ID number for panel data
	sort hh_hash 
	egen hh_index = group(hh_hash)
	xtset hh_index time

	local key /// 
		coffee_prod_dev_cat   
	foreach x in `key' { 		
		g l_`x' = l.`x'
		}
		
	* we now just keep 2017, having the lagged values for 2014
	keep if time==1 
	* keep l_coffee_prod_dev_cat coffee_prod_dev_cat
	label var l_coffee_prod_dev_cat "Coffee prod or dev 2014"
	label var coffee_prod_dev_cat "Coffee prod or dev 2017"
	label val l_coffee_prod_dev_cat `: value label coffee_prod_dev_cat'
	* transition matrix
	tab l_coffee_prod_dev_cat coffee_prod_dev_cat

*** Income sources as percent of total income	

	cap erase temp_inc.dta
	* reverse desired order due to append
	foreach s in other remit store labor ag_other animal fruit grains coffee { 
		use "$survey" if year==2017, clear 
		collapse rank_1=inc_srce_1_`s' rank_2=inc_srce_2_`s' rank_3=inc_srce_3_`s' /// 
			rank_4=inc_srce_4_`s' rank_5=inc_srce_5_`s' income=inc_total_`s'_dol
		g str source = "`s'" 
		cap append using temp_inc
		save temp_inc, replace
		}
	order source rank_1 rank_2 rank_3 rank_4 rank_5 income
	egen total_inc = sum(income)
	g p_income = income/total_inc
	list
	
*** Household water sources and use of irrigation

	use "$survey" if year==2017, clear
	sum water_source_* irrigation_coffee_prod_any irrigation_coffee_dev_any irrigation_grains_any ///
		irrigation_patio_any water_less_food_irr
		
	
*===============================================================================
*	Clean up 
*===============================================================================
		
	local temps temp_1417.dta temp_inc.dta temp_sever.dta temp_bean_season.dta /// 
		temp_tab4.dta months_temp.dta temp_a.dta temp_maize_season.dta mmm.dta temp_mo.dta
	
	foreach f in `temps' {
		erase `f'
		}
	
	
